{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import copy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.23699977  0.74029037  0.19926341 -0.44462407]\n",
      " [ 0.74029037  0.71792872  0.47324132  1.75475346]\n",
      " [ 0.19926341  0.47324132  2.06243397 -1.14174577]\n",
      " [-0.44462407  1.75475346 -1.14174577 -3.93667186]]\n",
      "\n",
      " Eigenvalues:\n",
      "[-4.83004377 -0.36887762  1.45552182  2.35009063]\n",
      "\n",
      " Eigenvectors:\n",
      "[[-0.13390826  0.93696931  0.28237085 -0.15628113]\n",
      " [ 0.32371695 -0.25620021  0.88539067 -0.21366365]\n",
      " [-0.17090423 -0.12251828 -0.20371314 -0.95618093]\n",
      " [-0.92090589 -0.20356628  0.30797894  0.12506833]]\n",
      "\n",
      " The 0-th eigenvalue is the maximal one\n",
      "Dominant eigenvalue = -4.830043768468989\n",
      "Dominant eigenvector = \n",
      "[-0.13390826  0.32371695 -0.17090423 -0.92090589]\n"
     ]
    }
   ],
   "source": [
    "# Randomly check the equivalence between the maximal eigenvalue problem and the optimization problem\n",
    "\n",
    "dim = 4\n",
    "M = np.random.randn(dim, dim)\n",
    "M = M + M.T  # 对称化以保证本征分解存在\n",
    "print(M)\n",
    "\n",
    "lm, u = np.linalg.eig(M)\n",
    "print('\\n Eigenvalues:')\n",
    "print(lm)\n",
    "\n",
    "print('\\n Eigenvectors:')\n",
    "print(u)\n",
    "\n",
    "n_max = np.argmax(abs(lm))\n",
    "lm_max = lm[n_max]\n",
    "v_max = u[:, n_max]\n",
    "print('\\n The ' + str(n_max) + '-th eigenvalue is the maximal one')\n",
    "print('Dominant eigenvalue = ' + str(lm_max))\n",
    "print('Dominant eigenvector = ')\n",
    "print(v_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f from the product of M and v_max = 4.830043768468987\n",
      "The largest (absolute value) eigenvalue = -4.830043768468989\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD5CAYAAAAOXX+6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO19d7wlRZn2U+emicQZkggDiEgGGREliKiICK6yiLqCsgbQxawrou6u7qfumnMAAUVFQQEVWBGQoIAgzBCGnGcIwwRgcrhz7zn1/dFdfarffit0OmFuPb/fzLnndHVVdXX1208971tVQkqJgICAgIDeRaPbFQgICAgIsCMY6oCAgIAeRzDUAQEBAT2OYKgDAgICehzBUAcEBAT0OIKhDggICOhxDPokEkLMB7AKQBPAuJRyti39jBkz5KxZs0pXLiAgIGCiYO7cuc9KKWdyx7wMdYxXSymf9Uk4a9YszJkzJ0fWAQEBARMbQogFpmNB+ggICAjocfgaagngKiHEXCHEKVwCIcQpQog5Qog5S5cura6GAQEBARMcvob6YCnlSwG8AcBpQojDaAIp5VlSytlSytkzZ7IyS0BAQEBAAXgZainlwvhzCYDfAziwzkoFBAQEBLThNNRCiKlCiOnqbwBHArin7ooFBAQEBETwifrYGsDvhRAq/a+llH+utVYBAQEBAQmchlpK+RiAfTtQl4CAgIAABnniqDuCt515c+a3Y/bZFie9YhbWbWji5J/dmjl+/AHb462zX4jn12zAB381N3P8xIN2xLH7boeFy9fh4xfemTn+/kN3xmv32BqPLl2Nz15yd+b4h4/YFYfsOgP3LlyB/77svszxTx+1Gw7YcQvMXfA8vvbnBzPH//PYPbDndpvixoefxfevfThz/CvH7Y1dZk7DX+5bjJ/e8Fjm+Lffth+222wyLrtrIX51SzbU8scnHoAtpg7jd3OexEVzn8oc//m/HojJwwP45c3zcfm8ZzLHLzz1FQCAs/72KK65f0nq2KShAZz3nsgl8b1rHsZNj6RD6TefMoyfnHQAAOCrf34Aty9Yljq+7aaT8J237w8A+OJl9+K+hStTx3eeORX/c9w+AIAzLpmHx5auSR3fY7tN8F/H7gkA+NgFd+CZFetTx1+64+Y4/aiXAAA+8Mu5WLZ2Q+r4wS+agY+8ZlcAwLvPvRXrx5qp46/ZfSucctguAELfC32vfN9T11M1Qhx1QEBAQI9D1LHDy+zZs2WYmRgQEBDgDyHEXNPyHIFRBwQEBPQ4gqEOCAgI6HEEQx0QEBDQ4wiGOiAgIKDHEQx1QEBAQI8jGOqAgICAHkcw1AEBAQE9jmCoAwICAnocwVAHBAQE9DiCoQ4ICAjocQRD3W9YeCewKLt4T+2QEnjkmugzICCgowiGut9w1quAnxzS+XLnnAv86jjg7os6X3ZAwARHMNQBflg2P/pcmV3KMiAgoF4EQx2QE6LbFQgImHAIhjogICCgxxEMdYAnghMxIKBbCIY6IB9EkD4COoTxDcCaZ93pJgCCoQ4ICOhNXPxe4Ou7dLsWPYFgqAP8EOKnAzqN+y/tdg16BsFQ9xrm/Ra45cfdroUFQfoICOg0BrtdgQCCS94ffR70we7WIyCgVyDlhPeNBEYdkA8T/IEJ6AKC7BYMda2Y91tg/cpu1yIgoM8RDHUw1HVhyf2RjPGHjUTCCKwGWLUYWHxvt2sRMAERDHVdaDWjz+cf63y5t50TxaBWCmWoa5I+lj4I/O3r9eRdBf53B+CbLwZ+/Mpu12TiIZCE4EysDUOTo8+xdZ0t987zgf/7BLByIbDTYcDOr6o2/7o06nOOBNYvBw46DRieUk8ZZbB+RbdrMIERDHVg1HWhEb8Dx9d3tlxlUG74BvCLNwHPPlJNvnWzGvVCE6FLBhAERh0MdX2IO5eNUY+PAuuWVVws6dSjVTsza2LUslVPvgEbAYKhDoa6LkgPQ/3LtwBfnVVvPfqFoSaGOjyUAQSBUfsbaiHEgBDiDiHE5XVWaONB3Lmao+YkC26qvxqVaco1PyzKUIeHsnfw9FzgC5sCj/21yxUJfSIP3foogPvrqshGh64ZnJrLrW3CiySfAV3H/Bujz0eu7m49wsvbz1ALIbYH8EYAZ9dbnYDK0S/Sh0J4KHsHRe7F3POAZQuqrkjF+fUffJ/i7wD4NACjx0cIcYoQYo4QYs7SpUsrqVxfo2cMTkUMuGPX0yvtFtCGZx8aHwUu+wjw3X2q7S898yx1D05DLYQ4BsASKeVcWzop5VlSytlSytkzZ86srIJdwdKHKmAFPdK5Kpcqal7rIzyU/Qs1yQuoOJop9AkfRn0wgDcJIeYDuADAEUKIX9Vaq27jhy+LWEFVuPM3nZt1Rw1dv0kf4aHsIeS8F3qIpW60S1cj9AnnUyylPENKub2UchaAtwO4Vkp5Yu0163fonesPHwCu/VKXKtInUR9JMeGhLIx7Lo6iNFY8VW2+3qMyafg7oCz6jW71EXqko1YtfdS+zGmPtFs/4o7zo88lXQrO0hl1pS/c0CdyrfUhpbwewPW11GRjQ8+E5/WZMzEw6hKoeeEsZ/F6rEFwJlaJwKhrQ690rqrrETYOCDBAN6h1MuqVzwA3fmdCGfBgqOtCr3SiXqlHwMYPWZNGTfvw704G/vJf0dK4EwTBUNeGLhnIjGGuqh5B+uh5qLarzH9cIuqjzvuoFhqTJSNLVi+JnK9zf166SnUjGOq60CsGp+p6dNuZ2BwHLjkVePbhmuvRj6hLo+5y1Eddz5La1OPOX9eTf4XoL0P90JXA13bu/GL8hdAjhrpn6uEJ10O58HZg3gW9ucXZmue6XYMI3dqAOER91Ib+MtRXfR5Y+xyw/Ilu18SNjY1RV5XP/JuAlm3t6R5pt7y45xLg6zsDT97a7ZpUiBLSRz8w6j5CfxnqvkKOzlUr++ihTv7Y9cDPjwZu+o45Ta8+lK4ZnvNviD4Xzau/LibU1Xa+DD0w6trQX4a6Vx9iDnnqWud1VZZ3BfmsWhR9Widk9Og9HpoafQ5Otqfrah/tdhx1n2nUfYT+MtQJ+iGWN0/nqtNQV7TFVRJRUKLt1T6SrXF3OT0H13rZPdQnq9Koeybqo1f7ROfQp4a6Yqx4Grj9F9Xm2S1GXVt4nkLNhrpXH8p+2IGmtroVkD569T5y6OV7GqPPDHVNDfrLtwCXfrj6jWa90Q/SRwXoZ0YtXYy6w7jvj8AT/zAc7Ba7r2lmoimvXu0rNSDXWh8bLVYvjj67NVyrtcP1kEY9MBR9lmHU3V5DpVeMw2/fFX1+YUWNheSVPupqm5rWr0my6yHZyoA+Y9R1oYYOlkuirkhHZvOuOL9SGvVA9FkJo+7ww9UPu6RX4UcoVX5NGnXdL8deefla0J+GuuqOmDjLq8y3W87EHg7PU9JHc8ySqIfqq0P2GKO2oqp+nDOfuqI+jHmVLaP3mbRCfxrqylFDWNPGFp5XRT6JRm1Zo6FnDaFDo+6D4XN+9EjUh1GjLjsS7dW+lkV/Geq6HuK8jqJlC4C1z7syzVOBHGnzooeiPlQ7WzXqkhhdDZx3bPVrgfRD1Edd/ajIhJdOMOqqJMM+eMn2l6FOUHXD5hzWfncf4Lv7ObLsVnge/V6j/p0Xqi51huctuAl4/G/AFZ8ulw9Fr0V9cJhoGnXVcwR6GH1qqCtu2CI3atTlbe8RRt1LMxMTQ23RqMvWd9pW0eey+eXyycCzXj3x0Fc14aXMCZ1g1BNHo+7P8LzaJZAO5yVbwKrFwFLHXneF6ld1W5XJT0kftnWEXfn7ylPz/dL5gLa7lAxr7YWHvocmvPQTo86LVYuA8VFg8x07VmSfMWolUVR9g+rIl+lc1/w33+mkBG47G/j12x1ZFuj8VTsTy+TnI3345m8a3idacoX3MpEUGunv6UTVlVcaXapL32nUBev4zd0i+bOD6DNDrVCX9FEzC7jhm8CapVxiYHw9MO5aZ9unfnWH55HZZ6Orc5wan2sNzyuJOmyUMgiuFfSqQJk1rWXVhKPEhJdOMuoVT0U7tcz7XXVl9hj601DXNeTp5sxEn2vySsMM06uEnt/NPwD+5wXtVfF8z60zPK+WvkEYNXtvK5A+FtwcrWl936Xl8unWrj61xVGbyovv9eL7os95F+TMoBfkKj/0l6GubdJBDdJH3jpaddsceWauoUZn4p2/iT5XL/HMogNRH7XMMvWRPirAwtujzwV/L5hBlyNTOrVnYqa8iuZBbFjj2NSie+gvQ52gD6QPW163nRMN1fTyvV4SBQx1nYx6bE30OTzV89wORH3UwagT6WNA/WBLXEE5BR9Ljsg8eRtw6Uc65N+oiVG7pI/CYYlavutXAF/ZDrj+K7mr1wn0p6GuzZnYISYy92fZ8n12VC4ij9TVVgCwYW306W1YlPRRI6PuhPTB9ZMqYpcTQ10yL70NfvEm4PbzgLG15fLMW26dMqJqnyoZtVo5864Li+dRI/rUUNdlUDsUnpfx97V6X/pIBh26oV6Tr4yEUZfQqMseLwIqfdQlLaj2UYtXFc+I+alAnfO+MDq1Z6LJaZr7BScMf/ce+tRQV8yaKveWA7k66jd2Bda5pqR75lm39KHXQUkfvmVUqlE7wvOqRGJAa9aoS0eXyNRHnFmJ+uSN+ujwDi+lfVa9FFJpR38a6rpCzjo14YV7dhbc7JFnkciQGp2JxjJN6Tykj9L3oJb4vOijU4y6qHGtbZp7l6M+nBsHxJ8P/Rm45NTy+fYg+sxQ16Qlu9YaXjY/2lEjX6aGnyV/yEujLsKo3ad4wcZevA01YdQP/B9w4Uk0UaHq5a5Lrjw7FPWRkVhK5pP+sVyeXuXqjLrSjO3l6debK0TPUclzXg9c9N4c+dWHMIXcJ9+fHAqMrsy3m4axjo5OZ880nT+nyWXKrVH6MJZpOpXMGrzgX4rn5SqjSuSK+ihTTllDbZHwCrVLCemjk1EfRcti20SPmLkFeBLA8ecUy79C9BmjVqgp5MzUmUdXFsnU8LOhjLyyhrHzEmZeZ3he8ltO6aNO1FJGDkZdxRT70jMgtTpkIiQKwNtJR4hEZcjBqHNly7RTj6I/DXWdIWd1Z2nqVDTQ3rQmiKuA2jXqEoa60BT4nKhjGrNX1IflQV+1CHj6do9yaoijpseK5KfjC5sC//dJQ/puMerCGfuX1WU4e4QQYpIQ4lYhxF1CiHuFEF/sRMWs6NbqeblmLeVl1DmZsG/n7QSjzhueZ8vTWV9Xu9QQeaDyKRr18YOXAT99tfn4sgXAM3dVwKhtL9MK+8FtZxuK7xKjLix96OeRF+3Y+mJ51gQfjXoUwBFSytVCiCEANwohrpBS3lJz3bKoJYxOz9+Rr2zCexBSWqN2MVdfQ92B0UdeZ6L1d09DbByq6saihWoGjXm0Y6b+LulMrcR2yCficrwrZqhCmVEPhx5f66Ow9GFpk0JyZ31w9jwZQS2RNhT/6/L4oK7iXYzaIzLDlZds8cd8mLCkRsiUv0c98oIOq33qYsqDQl9Nr/SElxqG3xmmy+RbicZZVdSHw0nmn1HxcjsR6lrnMqfrcwQOdABeUR9CiAEAcwG8CMAPpZT/YNKcAuAUANhhhx2qrGMW3Vo9zyeEzpWXUaOmeTskhm5JH5DAxe8HNp9lLtN4qiFdau2PstJHgReIC50Kz1N9oKxGnWqjCpyJ3uV3ej1qi9TjlS03Aog/+9FQSymbAPYTQmwG4PdCiL2klPeQNGcBOAsAZs+eXZeIrAqrJ/tOMGpjp/PQqItIH1U7E6UE7v5tsSJMxqIuRl31usydmvBSOjyPqd/qJcDYOmCLnQrmbciXK98rbYlyaSRL6ZmJMpuHMtSi7HT+apCrR0gplwO4HsBRtdTGvyY1ZeujUfvmlZNRZxPaz/XNv84p5MlPJaM+UjMVq3QmVsWoSRx17RNealjr48xDge85NmR2Zut6Pupi7RJ48tb2ImCZ8ooyauOXtqEemV4s74rhE/UxM2bSEEJMBvBaAA/UXTEruiV9VBL14Rsh4Yqu6DSjVtnV4ExMMWrP8n2ctbWF59WEquKo9XauMj44l6GusN+tWgyc8zrgDx/MVx8X9Cn7tF8pZ+LwNMv5nXPV+fSIbQFcJ4SYB+A2AFdLKS+vt1oOuBqoOQ78+GDgoavyZuw4nEf6MOVRgn36sMW6NGqrHuj78jE4myrVqDsgfdTGqMsanpqkQWXsXdJfVc7ER69Lr9c+uir6XDSPL68K6YP2KzXKUyGZrSZw5eeAlc9ky+8AnBq1lHIegP07UBd/uG7MmqXA4nuASz8EfOqhHPk6Gj6PRp03PM/nfC/po25G7VOm6VyDo09vV18N2rW5bZ56uUBXz7MuUFWmvZWhLVrvmkdPTkbtMeLzgU+ctk99XLA9T/T7/BujreeWamJCq1nBkrR+6M+Zid6dIO96ulUyal9Dmuf8AtJHT8VR6yxaa8s8GrWvNOKV2BNeUR8VbhxQ9J4lBLFqxu8ZOdLxZU6rYtRMHsnohJSVkukqGGF7or8MNW00Z8KiBRhQBaMuo1F7rfVRV3ieZVhNy7zpe9Hw1ZZOGgx16ZmJdUQeVBT14TtaqDM+mMNF7wXuvsgjW5f0UcNLMl0A+VrVokxCy8Mm8ZnOrx/9ZagVfEO4cu9Q0QFGzYUC+Z7v8yDU3nlc9QJw9X8Av3wzk84gd6RegDkMsat+XYn6sMkinoy0MieZJ+65CLiYW84zr9RAXpJj64AH/5yvLmy2NTFq3e/iLasYRoU1oz8NdbekjzxRH2VnU9XBqNc+H4U5lUGp6cn6g1wXo65xrY+yMxNrN9TcqKfCGZN5Z4X++TPAb94GLLwjX3FG/4NJR65iBE2uUX1uWA2seEqL3a6BCHigzwy1r7OlxI3Tb8Szj7Q9zkBkXMrEQQPlz3flk/ld+/6zo6MwpyKwPRRFwvOMjNq3Hh5lVP0g2TRq9duVn402RHDVzXa8sEZd1nA58s2rUT/3aPT3+pLrZiSGu05nosG3s+554Nt7gn3hBUOt4dafRqE6eWavKRSRPvS8f3AA8Kt/bn9v5TDUtTDqIs5ELd3S+/3KtuXrjO/2yCPzdw7poytTyHPGN1/1eXs+JrTKMmpVTsUjieYGYNHdHuF5pH8WlSAz+TpCUat0JppeChwmrPTx3KPAr45Pz0C69v9FnylmmydMKA8kMjfoSW1ZE9nMHrfmled3j3RFpI+qp5D7aNTGLAz6Xh7poxuTLnwe4pQxMhgm14Pd3BBnX6EzsYyNVPW48VvATw4Blj3ulz6pgqpHmd3B0W63zGCxKmeilofv6ME3TUXoLUP95zOAR64GHv9b+zd2uOnLugpsd28zFN1m1IXC80oYVt9zikgfJo3amzH7LnNaASijLuy8ctRnfL1fOleMftUTfRTWLHUkr2lmosmJX5kzkcmDfme3vZuojJqbzJA8HDm0x1JRHzaDlIdRm/JI/itwrnZecwPv3PSJoy7yINseiio16kqdiRUbrEbJPRNd9fFl1C7/RNXSh0LemYlFn0OaPimXvKTLMmqW+JhGjsGZ2AanBaqQqNa4f0cs4zW3MmoH405lVVL6cBnEb+8JXHSyPY2pvEIPcsXSh6wpPK8TUR+ufH1mTXIYH/VL55TVclx3npGSk0FSw1eR9NEVRu1hQyasRp1MLNBulGIxedaDSAWy5ym+Zc+7Co1atjzr5eG0u++Phvz17yUYMJdPVeF5KRZdJaOuU/qwMWqf8DxH3RWjzuW0K3DcOy1l1DmjPqpyJprKrTXqw6MNg/TBMOoq1yy2nZfRprS5/JVo1GUYtU/n8WHUZQx1CSnFJEvYNOo7zo+iflYtTtcjbxmlkJNRm4y268H21ajzyD+uF4itrAzL7NDMRGrYVbmdiKPO5E3rZGHhNaI3DTUYjbpIZEDuFzmjUTe0datks4K3uIO1p+rCnes6zdHhTHn71idlVNUhT7kho1GrVdks9/b2X0Sfzz+aroexjBoeHp8JL175uKQPpVEXvMYihiuPoebuvSl9yt9TllG7pI+ycefa3642NEUu1YweM9T0gUB7xbLmGJLG8w7RKhlHDaRXx2p1MDyvqIGtklE//jdg+RPp+rCd01NuoBo1u3wmrS952PNII1VLH767kBfVqJsVadS5pA9bWdRQq/vkc30VSh+SOhOJHSg7kzPzN5MnO5IMhrr9m+5MbCd05FP0DcuwXV36yDMz0RqeV1CjLiJ9lNGUzzsW+P4B6XOcjNpmqAtEfdCHnYsM8q3Lsw+nJTSKsfXAyoVcpnGZnWLUdc0T4PKyGBuTg82r7St0JlJyQGW4SpyJ5A+fiJcQ9aEbak768GTUud/kLkbtK1vAnM4VAphKR38rEoRf0OArKAeXVfrwHA5m4qiF9ndygJ4Ufxru5VNzgQ1r+Lrof69aDPxgNnDFp831+8WbgG/tztQ7r0ZtQO2MWh3WjrueAd8XK6Dde48Y9tSzVPLFQjVqyqTzGsxl84GHrsyOANKFpr9mQgQxkaUPSxx1LmdiGc3KYqirYNSlDH0RRl11HDXHLHzlBr2TtzTpY5xNkgK3KM6aZ4Gzj0hv0WRi1OuXR5+P32Cunj4LNVVtZagti8T7zEz0Ds8rGPVhc/ga88ohfRRl1FVGZ+j5FWXUP3ol8OsTkLo+KpeaXlIhjhpghx2dDM/jGDWN+qgkPM/ndI4J+5xXc3ieawhoHUobGLVN1srUX/uu9rVbeKehjBzRD9bQs4oYdVVTyI3H63YmFtSo8/Y344QXUk5RRj0Wj8DY0ZdB+mAJyoRl1MyN1cPzfG98LzNq74e8KumjKtikjwo1alO5nDNRGdfUdkgWB5HxNwBL7svmm5zC+E5s8F2mk6L0hBfPcnzTZl768X3qOKMmskNZQ93OOP5gnnuaJ+fbyLXZdTn0mKFmGr5RJDyvoLeZjfoYRJv5dZlRF5I+qmLUns5Eq0atMxhD1IfLmQiSB0BCKA0vDU460THvwmy+NJ+6nYmJFlvQmVh1eF5Gp3Vo1DSOuiijdjoTKzLUnE0xtWGSNkgffMOXcSYWkT4yUR9aE8lWeUZdxtD3woQXZ3gew3i5Y0a5xChSM/krhqePelzShyH/p27T8qUvow45E5N0rvzzMOoqJ7x0S6M2OBO5cvNAsWShb8VlyLNF6qDXqwPoA0Otz0z0ZAxlwvO4qA+d+ZWRLgCetZvScfXLex6bTxFDU0D6aI2Z06UmvNimkJseIKExat1QG14aLkZt3GxXOyfPWsUcvA110fA8D8NoMr4+9UhejAazURmjJsiQA4MEkhdJRBPcz42r39eMHjXUWiM1mDhq36iPKvZMpDMTS01YsfyeTej5G03SBUZtYsdU1zOl83ImMi9pdV5q1GMw1HlYasZQEy181aJoWvsThigRYxGe7e4M+yrBuPMY6rzSR2ZkReUET2SmkKs6kvyqMtRSup8bzlBP3PA8dQO0BmDD82pyJhqjPqpk1GWiPhxlP3lb2ilmzKeMoXbEUet564wlOtj+Uw/P81mUiRt6NuO6mJyJLNP3uC8mBqf6olov/dYzDfUuuNYHVxf2OHMNj9/Ah5BlyEoeRk3S2iYL2fIqa9Dq0qibXCSZ4eWSjA613684HZj3u2Jl50RvGWpu6JZo1I5OksqmoEZtivpIGZRuatSODvnE3z3Lq1H6sEoIWrpfv7W9CJHPFHLO0Ko+YdKoOcejcaSjX4Oh3sku5Fz/KjLaKZlO4e6LgPOOaU/3tzLqHPquyVj5atSS+b0ITM5dbgSeB4pICMGMNOi1M/1n6f3AJe8rVnZO9JahtjkTizDq3NIHo1FT51UVGrXX6QUMrO80V1v7tZrAkvvNQ2Q2ntSgUVNGbWRcFlkrM8zVGXWcf8PDmZiHpS5bADz/WPZY4kxktNpCjt6C6WhZ9/7e//wy0keTkZqMeevSR8moD3rvSzPqOP9U/3S0i2tBqprR+4Y60ahzrFlcKq6S3jB9qnMrezx3ESXYheu6WUOdk5lf9xXgRwdFxprLx1VGSqOmhtrDkZeHUSeGw+RMLCh9nHsk8L39s3Wguw0JA6MuGkfN1YVPkP769Fz/cmQLmPMz4Hcnu8syDv99GXVBQ+3a4aWsoVZ+J0X+UnVlfCFAMNQpcDdWaDMTfW98GenD5u2tYj3qMtKJ87odoXO2vBWeujX6XEUXJ7JIH6bwvKZF+tDh4yjmpIs8jLrMYvzUmZik1fpXNxl1Rju26eUSuPxjbRbuc90KzmVOqexUFaOueK0PdR8LMWqSbvp2+couiB4z1J7SR5Uzs9InZvNujac16tJ6m+/5RaSPcWQ7eUFnommI7JJXWjZGbSjXpCvrX2nfEELTTBvMCeDZtfH96cEsrSv4eTgTfZ1qeRm1HpnEHLbmXUj6qJlRZ/I1yXAlGfW4Hp5H28WgUVNssm30ufxJYLVj898S6E1DrTcKt9ZHbeF5jEbdasE4M/GGb6XXmUjlVYNG7SN9NOjCQQXyAcwdN49GnXEAm6QPG6O2SR9x/rk06iIjHYP00SuMmhrqqsLz8joTU+WWYdQExqiPotIKw6gzTWYYWdM2mToz+vzOXsA3XpSvHjnQW4aaDY1Rwv84vG98lVEfKSNCpI9rvgic9ar2Gg0KK55uL/zC1s3HWBcwsLLJsKu8zNywepg16sOgC2fiqD0MtSmOOqMhQvPaF3Amjq0HnpmXPp4xeKQOgvhLcmvUVUV9UENNXs5W40uH9J6LaAHa/dQmDy26h0/P6r6eMG3FVTaOesPamMyoUTonfRhsjFH2Kbkpgid6y1CzDa8MRJHV8xxYtQgYXU3OtUgfprU+aB7nHAnc9F1T5fzqRuvxzF3uc1tNZinOotKHYYicZ8ILy4SYjk0dxeOjwNKHzGkUWEZNWR2pozp++ceBMw+N+oA6PjCcLUM/xyZ9FGnToulcjNpqGHMw6gxpIYx6zrnATw4GHrueKbcTjDqnof7+S4G5P4PXPqxGPwBN15lJL05DLYR4oRDiOiHE/UKIe4UQH62tNlzDJwaiho0DvrkbcNbh5DzLjeCkEZoGANYti/6xdfNlF1q6R68FzjwMuO1s+yktT0ZdyFDbGLVJ+mCciQND2fMpo777d0zSSeMAACAASURBVJEBGF2FzMOus7Q8jJpGDiin6eiqdn4Npm76OZlFmbolfZDjmXA5C7vPo1FnjBWZmbgoHpGoUMbUyKpKjboCZ6KUwKpnotGuzZloqrNr38aa4cOoxwF8Ukq5O4CDAJwmhNijltpwjaT+bnpEBtBzfPDcw+l82SGPg1Gb3vhs3Qow6ucejT71YSaH1nh7WGcrz8uoGJhXqQkvkpcXaPutWx49RBvWMn1Cq5fK32etD2rouWsYMEkf6sVPJrzkdSbKVvQCLxJmmc7I/t0ayZFHozYM/03rptQV9VHFVlwqj9a4QaN2tKlJo+4VQy2lfEZKeXv89yoA9wN4QS214d6QSQMXmUJewcYBLS3Sw7QetWnmlKkM37r4/Pbsw8AtP2nXw8exVGiYrhh1ntXzmHbhWCudQq6+px4kz/A8Y9SHKzzPR/oo6Ux89mHgq7PcI6O80kcZllxE+siA8WtUyqjJSzpjJ3yIh2aoaRw1V0dfjbpD633k0qiFELMA7A8gsxqNEOIUIcQcIcScpUsLhqlwkQUJo9biqKvSqNnymTdpSp8ty6gtx3RD8cQtpF4GnP1a4M+nRyMOnS3YzrXVgeqwyTnxZy6NmpM+OEZN0rV0Q214OIUo6Ew0sU+ZlmUEY/yT5QzivBbdDXx771jm0uUG8Hg21t0fvsqQgNbJmIB8NbxU2cqUYN+umYmZqA/1Z05jVscOLzqjtjoT+a/GOOpeYdQKQohpAC4G8DEp5Up6XEp5lpRytpRy9syZMwtWxyJ9pBi1KxtPjZo7j2XUis1VwKi5STUKAyPtv//vE9qmrZYLVnsBylZkQHzC87wYiMmZmCPqg1vX2SV9pBg1N8lJZ9Qu6YPpRzbpQ2f7utFW5zYIo154O7DiCWDBzX6MmrbpI38BfnmcO13muItR24yvgznOPS+6Hu5YZmaipR466alKo6bl5NKoNcJXSProrkZtEOXSEEIMITLS50spL6mtNqwzMW6wZh5nIuPs8asAWEatD/u5hyCzQH5RRk1kgfFRYHiqdq7tAWzGbGEAmHUoMP8G8zm+zpf0D9FH2UWZOOmDOhNVe6YYNflMjpPyMzqpKsNwzbreqo9o9HoapQ9ovxcw1BecCIyvc6fLJiBfHUbGlpaWddlHos8vrMieSzVqhXXPAw9fnW37qqQPuiBSxrHs05856SNPeJ5B9ukVQy2EEADOAXC/lPJbtdaGYz0Ju9JjlWuUPmydXrb4sjPLdNrqJ81Mf3CEJKXtYclXMX8xAJz0h+jvL21lqG8Bjdo64cWTURulDwOj1h8OGrWhH3fKHfrfJkZNpI9UPamhJm0gGiRbzzjqgeFihrpORm3Lh4ujBoBr/jv63Pcd5Fzm5eoFUxw1zbeo9OHBqLkwXfr78LSOGWof6eNgACcBOEIIcWf87+haasMy6vjvDWvhPZRKpI/cFbAfNjLqPAtG2aQP4sxKOpKns0SF5w0Mto0+y6jd2WVZm67T06QG48h56znpg27F1bJIHylGPcaUY5I+TENXLX/dUHP7MCYTXsgLSDSKjVK4UEVbXdsJ0nnajG+p9agNDjUT0dB37amUUWsvWc5ZnYd4pKI+LHMzMtdOWP3hZwAvPLBjzkQno5ZS3ohOTb+xGeqxtfzDasunSPkuecHJqD1fIhwoo/ZdfQ6IOozOFgDzkLwIo/aVPqzOROmWPjJRHzaN2sWoGaOdkcV0Qz2cTa+fw647o7LSyjLOTCQPtTHKJAej5gy1VfrIwahpPpRRU6j+1xo3MF8J3PIjYO8TgGk5/FjG+1vAUDfHeGeiqw0TRq3Sicjg9xCj7hw41pYy1DkZdRVTyHV4Meqc+qKODKPWlmF0nStb0b8UYxV8fax1NEV9WAy1icWyGrUr6kNn1BscjNqhUVulDzosl7wDUf+b2xYOiBl1AY3ayKhzPPzqvqd+89SoWSNvyce11geNOkqaNi5j8b3AlZ8FLn6vuUy2HobnKw+j9o2jpnVOzidRH0L4j6QqQI8ZavXJ3IyxdTyrYvNRQ9UK4qhTh2tm1NRQJ2uIeBgBxaj10Clu5wrf/EzMy+lMtBhqyOyEHHqObjxSca7MQ8kyatNLw+G1p45O7j4ZGTUdueTQqH3SZY7rZTmMrS1vzsjb6kFnJnJyQWqndvKSVYZxPeOoTGeU/poYZIAdsXm9JG1x1AxBM2rUmm0RDeQOPSyIHjPUNo16jfbWq1H6sDJqgzTSYuprLMOSv0n68LluXaNOIMC/WDzax7TGsYuhWzXqFtgYXBujtkkfLUajtoXfccf1F4BuOOnLA9DWfGEMdR5GrdK6JtiYM0qnzcOoXeea0gIaozYkb45p0ptNo/bw46S+G9aPyeNM1IkGdSam/EaG0SsnfTQG/O57BegtQ83GUce/ja3jjyssvANY8Hdy3MKoWSdbjzFq00IwHFiN2sCo8wwVk3M8dU/XhBfOUBtnJjqciWPrs+UbnYmGa9bzdUkfyYQXyqIML8TUcaYORaUPF6NOfbes9eH0yRiMVXI9yB5PptkzjNoXJkceJFijnac/s3HUHs8ILUuI+PmakNKHhVHrzkSuYc86HPjZG7Lnm2D01ro0aibvXBq1pYwMo46lD9+hnQrPSyCieOrr/zeddsNa4ItbZPfa02HbmDZTdg6NGgJ4/3XmskZXAasXR3/rGmKLuffrno9/MrS/VaMmxpMyas74m6QPOhFKMe8VTwM3fke7hmb6uIlR51nrgzW2vhp1C9ahu0n6MEmKzQ1aWKONUTskSduLh4vTz8Wom/xCVYvmkd8MjlTdsSwavRP10VGwhlpp1GuBwUnZ42w+ZKhqK4ue52LULEPNw6jzOBNJ1IfrJUKlDyGi/fSengsc/pn278sXRHW+9kvAnm8x5MfoyyZwIVNcHlJGnZu+kPR0V57R/pudQq6dt/a59LFFd7dXw6P1yvQthqmnGDUjpySMmhhq09ICPz0CWL0oWweFuqI+vJ2JDo3aJH2YDO3oqii2eO1zBkbtKxOYNGpyrT5RH/dfBuzwirQzkY5knrkrXkZYr4LBmaiP1ids1Me0raNPjgnJVix/wM0wvaQPk9bqMIbc8TyMWrFKDhlGTaQPl/ygL4oOmMuhQ3kdiQ5Ljay56NwaNS3XNGvQNYV87fPtclpN4CeHAA9cnl3ljqtLUidNk3RKH2rIzLyAUvcmbkPdSKfKipFX+mg1gQevyLJ9G6O2xVE7nYkGVmkiQBtWazNptXrR++ty8mdePPG927AKeOLm9u8uRr1uOXDhicCvT9AY9Vgx4xqciRo+Fg8/jNpihc5ENk0VjNrHEWRIY4r68JE+1MzEFKM23N7EaFkemIyRtXXIHNKH6uA+ebOMWusDifQh0/VtMIba5EzUtUdf6YMyam4vTbrrTyZP5GfU/zgT+M3bgXsvSae1atS2vClDpU48yioda31sWNM21Gwctbla1nL17xeeqNWHGYEDkdQ3urr9Ylk2H6moD5+Xk1Gj1vpCI2bU+oto4R3mvEugtww1kI1NNBlUG/S3nitN6jeLEQXiG8Ixag+nlc9xozPR41y11ocgzkRTWuNxA6O27UBtZNTMOSyjthlqlS+ZQj66QjMAxFCKgXTa1PlJonR9ZSsbnkejbUwaNTfSYmdw0gkvDKNuDJrv88qno8/lT2p5EmM7NMVBNiyM2ldvNfWr0dVR+VHm6bZNZ2DP3/QyzaRrAqsWA/f9Mf379f8TbT6hl5c4E8ft7WOsAyd9NLJteP4J7rwLoLc0aiBrqNX2UnlYq8tYzr+pzcbSJ1bAqEto1MaZiR7ec06jdu2GzTJuNVzVjGyrVY2hNkofhryb4237t/i+ODwyzn+tdv+oU0y1wfInon+b7aDVy8AapcyuQyJl2rPPbbQMZPuFEPwogV4nx6gbQ+Y+ogw7XfeGGmrfmYnUyHhr3R7Sh7QZagdM0gdFaxz4/an8sTGyhoovo3aF59EJLy1D36sYfcCoJTAyLZ3Gm7UaOtTPj04PoVLndUKjNiDDqKn04XiJ0F3IdeaTYv0+0sd4+m9fQ+1aPY/Wi56jQ5c+5l0A3PTtdlsoR6I6PyV9xN36b18HvrN3nMYw/NaNSWYtb2LcVb2pRs31C+4+ZybKMO0/MGQ2TKp/6LKKMogHfww44ykmJNOmURucc7ZrsEKmNeo84XMK65YDSx9M/2byYbSa7TW+KZqj/EuoNQbrc0TTJ2V5MuqaDHVvMmpq+IankdlMFWjUHK76vCPfJl90HkZtq3sZZ6J6szdIeF5yXDe8OaUP2TQ/LLReVmei5Bm1TaPW837kmiiyAwDG4xjqkU2zD0tmg18w90W1qWZM2N2vB5GRPjKMuoWMceRePnQmHZemMWh+cSlphq5RIVtR3xmZzsySJMgwai0t3aTZmI+lPyrpg2XUHgby7Nemt8ejdVZQ7bRyIZ9Pcyzdr1T/dTJqVaZBfkkWpmq0DTXnH6kYPWioB7KdaXhKOo0va809hdyBKuKorYy67MxEqlGT4wquVdAAwoybOaQPjwkv1JAapQ9iEBfclE0zMOTHalzSByTzAiGGhl11Dfw95e4VF39NMTAc5Tf351E46r5vbx9T7ZRxVOp1d0zCsIX2feNF7msA7PkPq9Evw6h9wmapkQYs7WQY4QJRG3ETZZq+htrAqHWSo6I+JiyjTnWmZlYScGrU6njFhtokjTidn3paiw6+zd7p73StDy+N2sSotc7uWgUNYBh1NzRqTfowQUkFLlbj40w0GepM1AcjfRTRqLlRysAwMLYGuOyj0XfdUI+tjT710SX1NxjXd1GXopW5aiGwapFf2tTvFilOJ1U0zrloKBs3whgYsjs/mxsIadA1ah9nK2XUasKLJhuqqI+gUSP6O2OoS7DWMqh79bwtdwFO/lP7O50N5Yz6YCa8JHXUjESy16BnHLWLUafW34jbojHoH55neuCef7Qd6WBCYyiSQ35/Svs3m/SRCUHT21YAr/wwsNNh6Xq5oj44ZseOvDxmew4Mmds6MdTLs3km99o1nV07duZhwNX/4ZfW63eknYmmnVlyr2rJGeoRe59sbmiX12oCN8UzRH3jqH0ZNbUJE8dQk6EbZ6jLaNT/u6N/Xejayca1PkyMmumQnBaaJB9IH1MG1ccpo+/wwpWvn5tMXLDc/oym7RvyFf89MJI2wGuejddZyBGe9+i15jIVVCSEnpZdoY+E9yX11XRa0QCO/BLwYrIUgTpujPqgGjX4a/KVPkz+ABXJsE5j1Emopc6o+dOTuvqirPSRYdQFCZRp5JGZuauhuaF93vrl7Q2FTWSLwqlRK+lDwjmaqwA9aKh9GLWvoWbS6WzEBVqu6SabQgfZITjM9W80kDKuiUbdTH9y0PdMVNANYkpzdswwo+mbzOSNVNmMRq0zw+YY8PVdgMV3I/GW07oXBReL7OVMjKHf08TY0b0RpeF3PQ/ygmYZNRd/TWCK+ljxdDs2eNQmfTSAB/+khafRF1M2ayOM0oflHJszsei6GFx7DA63Hcocxjfw5zU9oz5M4Xn6ZDEx4aUP4q2lD6Ov9OHz5rQhtTWTeoDySB8GRm2Ci1HbhnqtliM8z1P64NIrvdQEbkQxMKx1bi2vxqC/Ru0DbscY1plocMjqjkjVXtQgu5zTrPTBGQkiUxgZNfP7+W/lGSRl1CuejOYIXPm5+LiHJGNEgaiPRPpg+kRhRm1oJyejNvgJCkV9qL6sjUaVL6IDzsTeM9QNJuqDhq15Sx+0sXN2FJ1RD454Mmo9TIx7sC11Fw5GbYNaPc804UXPI68z8ZG/OMpmwvMGNQ0x1ZEbTBx1CUPNbZZrdSZyhotEIzSooSaMmoJ1JjJ9bcOq9Pc8jHrts+ayo0LTvy+bH/+RCRzn8+FQyJkYG+qUM7qkM5F9oRnWSVFojvLnOTVqw0uIi4FPppBPWEZNpQ/KqD0NdWaIatjy3QTdUKuwqTonvDQIox4njNoGFfVhmkLORX1YpQ+L8Zy0Wfp7KtJDlz6YkYAYsERXFIA3ozb0iRSj9pQ+2LyJTs+RgnXLos9WE1hwM3/dJmfi1K34sk2zTE0yTd0atXpmaNSQnl/esFn2haaRt4P+LXu8OWZ+bnyIQeYaaVuYJrzUY1J7NDyPaJ5FNWrXW9EF/QUxOInRIkl59G9W+rAx6oH0Oc8/Fn36GGpuZqJpwotV+lBRH5YyJ22a1vrH10dTuqdsoUkfGqOmzpYqHS6sRm1xJnKLTdFoBLpBgB7JwoE6E+nDS/HoNdE/DqaFmqZtBSzmyjYZaiWvGKJcvFAg6kO1Eceoq9So9Xaa8eLs8XEDowbskonPLGCg7RRXTvz2Aft5BdGjjJp0+qLheVUy6sE4wD5PeF4R6UM/5+m58RZkvoya2eGFq1de6YNimEzpf+By4Gs7xeUwzsRUmzCMuhSYa2ClD8PwW49rTxj1AH+OyVBnomKk3z3jkFlGIG7DKVsaytZmyumolVH7GGrmmSgansch5T9i7ose9ZE55mEHXO0kdGeiLnvVY1J70FBXGZ5HGXVOLTQlfYwwD6QqJq8z0RT1QRh1awxYdI+diRz6yXYdbBq1Kepj6YPA43/L5msz1PpDklkfI762QQujrrQzM23pkj5M7NckfTgZddOcZ17QMsbVGuymYbyBUZuWIs2lURdwJrKMWqUvEd0DAC8+qv237rcyGWpTedwStBTOkYeKXpJE2ptIjJq+jU3Sx9rngecezeZhkj5yM2pd+qiIUVulj0bWto+vt3fwHV4Z16EF6zKnKelD81z/8EDgvGOz+VoNtXY/KHtNNOrhtm4+5xytTjkN9U6vsh/nDCIbnme4R9aoDxIHbJJs6AucxtYCwIzd+HMpTGuSj48CM3cHXnV6+rgpIiW1GziT3gsFNGql0bLOREf0jAt6v0n1QROjNkkfPnbAJX1oYaapZ2XCGOqBtKG1SR+/fRfw/ZemJzuk4jfLatS69DEJxvWojUuw5gzPo4xa5W3qcO+5EtgilhzUS8gUR52K+rBodAo2Fq+3i2nFOTUzce7PgL9+tX28wUx4sSET8UPAsjq9DQay6eiLlTLqhnbO5Z8ArvtyOi8A2HQHLQ/iu+AYtUl7pqDpVDx0c0NEFrjYfr3uClVo1M6oDwaqjarUqJO8dUOtSx/MC9QURw04+j+dRWmpi3op6XZlQjFqamgzDqP498X3RJ8rn9EO6Q+JIRbSF3q5AyMWRm1wJpbVqFXepk6jdpkA2p0vz1ofvnHUFIO6odYn1bTiDU5H2oZaRTok6XMyai6qIwWmPfX8Oc1UH1npbJiTPtRogM6onLKFVgWZ7heS0ahd4WSmdGpSx/ho1K4ZQ23QqNuVIV8r0KhNazYDbbJhY9RFURmj9iAq3tIHJiqjbmRvLGVV+qQKID1DSTfUtGOUivoYzg5xkzIrnPCSh1Gr2VGAZnwNq+d5h+fF1+ctfWgPSWs8MiiDmqGmyKtRc3HSqepy0gczkqDO1IRx6pq1kj7iT/0lQycj6W2QuT8y+6B7G2oqfcR9u7khatfMs+Bi1K4wMwtMM/9sDFm1k9WZWBA6CXE6Ey1RH0XWo6YQzLMHTCRnImOoM2tuEAakvyH1ELo8GjWna+qxmgmjZs7Nq1Ebp5APZM/JhP/oddYYtdIyfSa8cLtJ0/VEfJ2JetzoumXAmqWRsUnWVabaKbnGLcnSmpmyHJKBy1BzkSd009yMMzFu03Ne1z6HylKp/RUZZyI1YoWlD51R55E+KtCo9c0ZUnlY5IGOMWr92WRegs2x/BPcAK29PMPzgLRdqUn66O046iQsyrTzRtwoKUbdzHaMJ26JYi1tjHpgCBinD5fWPIpRs3HUORh1lIj/mWPUrfH0Q6+3j0DbqKjogFREhkn6YOKox9dHs8pM8cY6TBr1N+N41k22b6/FQDuuupfqOpIdQQxwSR+soebkI4NGrzNq6kxM5UkZtVYvzpmY0aiLSh+6Rs1JH3knvPhVAxvWRltrsbBEcSj5zmqoizoTdUZtcWgD9jhqH3D9avp2wOiq9gxTdjXFiSx9ZB4c8mCN0wdP06hbTeDc1wO/+meHRs008Ba7tP8enITca32wyoJNcxZZI0PXgh4g+rDqpOtXRp+pGGctr5t/0P6bkz6SHc99GLVBo1YYHI4W56H71unp1UNHN0vIlOWSPhwaNYDM+g7cDimA3VBTRp1xHtLwPGIknFp7DGqoxzTpY2DYX/owhef5slrTlHU9S6MjV6R3XqnKmagb5EGXRm2ZmegD7oU2MNjeFlBNIQdIeN5EkT4anPRhiNVVN0Jn1Pq6C1K2jcWiu+3GhxrIgz8KHP6Z9ne1/GSetT5MMxNNxkWI7DmUxVMjqTqGWkxeZ6j6NT3IrHOtl6XaKWlbX+mDeUgGRoChyfH6yRZGDbglAadk4GGoqc5vNNQk6oPmqbenWhta5Z+Z8FKUUVONOr4vSvsvG/XhS6nXLDUfc2nUrTHgyX9o6YlU0pHwvJKM2tSvhNZ/OWfixIz6MDBqdePVA2d68GSrbcQHhuzSBy1jtzem2d7giJlRj28AvrAp8PcfeGjUhhC/pANQRk30ztQDr73VRx2MWkcifeiMer1WP7Q73wm/zJ6fahfGkA6OWBg1NdQOA1ZE+qDn0F09mhuIs80Q9ZHKkxiDDWu0/BlGXVSjpi8J1WdNjPquXxvqbHAm+jLqNTGjpuu6RJmY8+Jecon/o8rwvBJRH1ZYri0V+aNFfeQNUigAp6EWQpwrhFgihLin9toAaUNtYgvJ/mdxA5miPiDbzGdgyDHhhTq9RNpRNmBZPe/pudHnnHM9HgQLo+bqQZ2JuhHS3/CJ9DE1fZwD3TkG0Bg1afvNZ2XPp4tVUShDPb4u2x55GbWLoZiGqF9YAbz+K9F32oamFztd6yNTb60uuqGmL3BO+vA11KYFqxJGTV5C918Wn2eY8JIJUfU0XopRT2MWg7LNNLTurkP6wpxzgR8e5G/ojIbaMBHJZwYihZTAU3MMddJsgmnCSxelj58DOMqVqDK4GPXQlDTLANI3JBX10WprfI0h+xTyjEEg39XMRM4yPD0n+tzhIGKEDYzaNkmD1agNjFpwjNogfejI7P8Gs0bN5ZGKL2cM0MBwJH3o+SZ1IoaaY+R5YFuwR4+hzhX14eFM1KUPmj8nb/lKH5wzGWjHp5s0fZMfhxpH37kEKjSRXWMkJ6NO+pR+D8aByz8OLL0/ctD5oGFyJpK2HZwUfer3yBetMeDs1/Cb7OrEyMSouyV9SCn/BuD5WkrnoE8hfzZurE1e0D4+NLltaDhDnYr6gMaohy1ebDCMhGqrQ2ZGrTB5cw/pQ8IqfbgYtf5Asow6h/Sh15WuKWGbSDHoCI0anNR+YdCHJS+jdk25Z2O1h9Jl0JmDmXBOH+mDMOpRrS9xzsSM9OFpqKm0oe88PshIHwqcAxXItp+voVZtxJVn1aiZtuOiPhSxMOXDQc/b5kxUhtq2A0wR6NKH7shv1R9HXVl4nhDiFACnAMAOO+xQIiONUT9xc/S5w0Ht44pR6yFR45oWSqUPdbNWLQQufm+eiqS/NgbAOol0ZHaPMEwoYZmINqTS8eQ/gMeu17LUO0JBRs09aGMGjZrreC5G3RhoM2qqUwtynb6SAAcxwBseFSmiLxBEnYmpCS8+jJo4E5NwSG6kxTkTHdf5qtOjENKMoW5GTLQ1FjNqg8E3Geqi0odiiVy9nXHUND0JtwXIKNjz5WEMzyNmbGhytAxvEUatY/djgc12BOZdGElBeh8QBkbd6+F5UsqzpJSzpZSzZ86cWTwjfa2PZ+4ENt8JmDqjfXxoctQwOivSjUGrmTY6vjeLBsdnlBAmFCeTBzHUJmdiHo36zvPT3+laHknUB2OoXZ1Gr0fGmageKk76cGjUsmU21HkZtXXtY4OhVow6YTyM9KHXlUYjcMbGpINO3jw70ioyhXzqTGDv47Mab2u83dcHhvylj+Y43898HXrJRsTcddscbh4atRDpfTh96rTHmy3SBykzkT5KMuqBYeD1X25v2qDPRjRNIZ8wUR+QwIKbgMdviBp60qbpw0r6SBlq4kxUEoeU/jcr01koo1Y3haSjs+B8wvPyRH1k0unSh2i/2cdi55YPo07qojOcPIxaN9SMAZISGDJIH+o61cNVWvrgNGplqDVGbWJz3FrCpuGr3p4vimctNoaQieThYuWLRn3IVtuoDY6YNX02dpyL+fdkr62x6NqMkUvgDTUb9cEsd5uXUb/158SZaAkRNZGEvEjIxGD7e0PrIxN6CvmBp0Sfzz8a3cDMTYilD71xUjGtrbYMIFtpWcQG2ukymrWBUeuODG9GbevgLkPNrI6nfhucbF6UiUNK+lgXRRCodvV2JnIMT1bIqC0QA8DOr87+rvqMvgGAlO26ppyJGtt2GWq9Pd/+a+D0+dGDyzFq+gJxhRkmQ2qGUasJXXmciaa9Ab016rH43pjIBsxx1DoGJ2dJAER27oMLutQA2KWPRKMuaajVtSc+D60OqaiPHpA+hBC/AXAzgN2EEE8JIfIIvfmx02HR5/ho/FZn3pYu6SNx9Ej3W/UNX48eukxnMUzUyBhqfVEiOsXcpFF7MOrkwTaFXTF1o9OxnYxaq+ulHwIuPBFY+kD0vRSjbpmdiXnjqCGBmS/how9EAzjm28CebyH1o9LHePRvKH6AmxuQGr5nXtIGmUNvz8HhSPZIpDoXo3Zdp7rv1OA226MlbsILVzd1nmsWrQ0rn46YpG3hLh9GPTwl6/+AzEZqAVn5kcJ3rY9E+qiKUWvO6ZT0EbeNHk3WLWeilPIdtZRsgmrk+y+Ldh+h+6ENTWGkD8KoE+mj5b5Zsw4Gtt4zK334MuqBQUC9UL0YtUn6IDe4MRhdpxBpw64/yJRR5zbUloc2mZrMMWr9IXFp1BVEfZz2j2hK8rd2J3k1IoO5xc6kfio8T9Oo1y2Pull8FQAAIABJREFUtMb1K+LdP3RD7RH1ER1gfmogu/xtgZmJ9F4qtMajuF4A2Hovf2di08SoPQ31fX8Epm3NH7PGUZN6DE/VIor0KCOdUStpxBFP7bt63lBFhprmr0d6GKM+JopGrbze828A1izJvqETRq1LH3rUR7Mdlyk9GLVp8RpfRq0/nz4a9Z3n86uS0agPZWwybI9o1HrdlHG0la/D66EtEEctNeljlIREUmPo27G5IX/CzkkfaRCNWjajtSs22Tb6rhsxNurD8ILlfm8MZNlrHpaeHDc4MltN4PG/AiObAtvtZzHUTPx1GY0aMGvUsEkfJP3wtGyMPgSvUbvq5jvhZbAijTrZ/1OVpUsfWn36KeqjMjQG0nqeS/qgU5WbGzQWJ92xlD7DXEDTo2jn1B4G+nCkvMQO0DhquhjRtvtl60uNHjWaLiOormWzHS31Kih9KGcincyQ0eJdHTtuT86JlqzNQdpqgGjUrfFoSvT07aLvzQ0GQ22J+jDVVQxkGTWnUTuHxAaN+pYfAff+Adhq9/jZMAyC2fA8T0M9ZUvg+HOzvw8MmettW65Xx/DU9vOZtAmVPgyG+hMPpL97h+dVpFEn28qZpA9Go54wjBqwb1w5OCmWPuLGGZ6WHl6v1wLpfaQP32GuiVHr0gynUZseLIU3/QDYZu8sm9M75dSZwOTNsvVNDJUpgsLRaVQ0wch0cxqnoTYwahWts3oxyU+xYEv+HEzx2kBW16WMenw0iqtNGPUocYj5SB82Rs1MeHE5pynovYx+BNY9HznH1TNhckqyUR+M9MHJFYOTgb3+Ofv7wBCMfYgLQeQwPDXrTFxwE3DV57W8DBNo1P1S0Nsm9Te1ERUxavWsm6SPhAjo9Z5IhloHvQmNobT0kSw/GkOtIjc0xVP6MDSsUaMmOtr2s9t/cxq1y1APT42uiTJNfSeUxlC7fC6qQxgMtcs4PPdI9Dk0xZyGy2PQYaghI1Yzdavsw0yvk8t/1yMjBy8AvOSN5nJMum7iTIyPr14SfSaMeqxdL98JLyYkE7Qk8JJjoqglbq0PZzQPw+b1vpMYalPoYPz7cWdHn0bpgzGupjyN0gfA6vAchmJGTQ37sse1OvlKH1pddJnPqFGXjKNmGbVWB9Xmz2vXUo+d7sGNA4B0B6DD0IHBtPQxY1dgxRPt48pQj2wSdRDX8Kcsoz7kE1EnfPRaXqM2DqNVEhULTeOo499kKz0EtTJqwrZsBmfrvYHFd0d/D9sMdUMzRjFc0ocqd7MXRn6G1DHmOile+eEo+ucLK7TzDLIDkG1j9eCqT8Xqp82MfktJH7rBMcgP7QKZnzRnYhKyxRgxei8GRtKTPriyGwNtYuDa5Fflv89bgYeuABbe6R+eZ+onNkb931vwv1MMTwUg021uqpPTUGttozs6TRNeykofVKPWn1PZbJf70BVaHSdKHDVAohyo/jicjvo4+KPp44mhno6MFsbBZ3IDYBjmIFpI/OivRbMnueGmupkn/DIbnRBlnB5S6fXSJQ2OPSYsbLCdToctomLSJu2/rYy6wRsZBc6A/NMPo89NX5g95sOoXaOQJJ1Jo1YMKC7rT5+KPqfMiNoko1F7Rn04nYki+sdq1ORcxfrocdOL3bXBgl7nxlA+jdp0vQ1TeF4OJGGa68zO6yQ8z9eZKKLQSBOSiCMPRr3d/sDnl/DHEukj7k8D2si31TK0zUSSPlKMmpE+ZAu48dvR9+Fp6Y6mDPWkTaIHhhrqmS8Bpm2j5efJnkyMWndoUY1aIM0gh5htp9RmAXShd10D0zsIN/RSBpmyW1tI2IhmqK3bYcX1M+VLyzj6G8DmsXNy0+2Z7EiXY8P/DC8YJYfQvFxRHwqTNkWyLrm+/oSPoZbKEBPozkTFqH0cbfTlyF2LnocvowaQbCzsK30YWaDhmvNA9a3xUXOb+DLqBkNSKLbcFdj7rdHfPstHiAFz2ybSR9yPhqdpjLrFt9vEYtQWQ60a7fG/xt+H0ml06WN8XXrdYAB450XAMd9qf/dm1MrDazLUAwbpY7D990wSE65+141yopciLWmYnE1Au6NlDLXl4dYdiC5G/Y4L2lOmaTnUqOr3govDzbxwHOF/OpRendTNIH3QCS8KQ5Oj/vGPn7TXW+aiPox7WhoYdTLhJX7hcho17WeZNmcYtd6Xchlqri/GyMOoTSGJeaBWcxxfZ3Y++uzTCWgvM0ud/vVPkfMdSEd8bbUH8OHbzXmy9SKMenhqO71sZgnCVnsAh33KnF8J9IGhpg8hHd5TQ708+py0adTQD1+ZTt8YJPKBQ6M+/Izonxej5pyJmmFSb/pUfQbS6fQIhJQzUbR/p2UnjDqH9KEbQxujFg1g19cCJ16k1dliqPV8OeNCwxDzMGqubvonrYNJu9ShzypUdWGZnYlRN9pRH8rZlDL+zEgIYBg1p49rL4w80sfAkJmd5mHUsoXKGPXYerNGrQy4U6ZUZMY2WhHt+68z6qkzgC13YfqKpa+p+qo0w1PTk6jouf96BT+KrAC9b6iTB18xDsoah9OGep0y1JuARWOQDKEcjPrwz0T/TBq1zVCnnIkC2O0N0ZT1VP1H0kY56YTauQPD2gPPDP+UAbIZTQB47Rfbf+sPct6oDyujdi2BqkKdmBeP7TwODQOjTtZm8DTUNOrDxOxYh6ZyJrba0gc0jZr2X4Vt9wFm7t6OREnqrrVHSvrIsRNOYzBePa8sozbpsDmgHNXjNo06rpNtn0Yge59f8aG0jAlE9eVeanqInQ5bu6p6paQPJVG20uce+aV2CG0N6E1DDcaZ+MkHgNNuzU4EsUkfHBSDVVA3Tt9xPDqQ/qqMGZ3AkdGoDRNeVJnUoTg4EpVF2QLVqNmoD2WoPTTqwUnAIR9rf9en4Nucdy5DzUXlJGVyDwz1CeSQPjJ1U4ba4Ezkwrbed036t1TMs2LUpqnMBsdn4rhTzsRWlo1RAzF1JnDaLdESBoB2303SB/OSedcf27/7atQ+077biVGeUcfSh41Rr18ZhVDSmHsKfS1oIFqC9FMP0kR8/6ELdSXPlmWkol4s6tyhKWnpQ78nm+9kr3tJ9Kih1qAaafo2wMzdsmxrcHJW+hiaYn7YTdLHB24E3n+d9jvpoNNi3Ws1DTfTtEWXRg1kvf2DI8COrwB2fGX0XWd3STSHwVAnsxiVRm1gt9vuF62XoUOf9mpbY4F7iOm+jaZjrPRBrqOU9GEIp2P1fEQPFt0DUp+sos6b+ZKI7dLIAjamfCTyg6x5rs2o9Zhh9YKn56qJUkn7SabOuvTBtMmmL+T7hQpB1F/GCkU1atfqfyYkzkRm/0yFP/4b8I1dgVXP2PPymeWrwl1Nz0KDvNxtpIDuzTo0SRtZt9KG2uVDKIk+MNSGYa3C4EhW+hieBiMT0CMogHYnHZ7CsxYFtXj46kWGeg5GCwfprECf8JKwX7Iex+AIcMTngdf+V/Rdd2zp0gc3zFd1N81aU2l3PyZroHRDbd1c1MF4TfowYBiCekwhLyt9JKyHWRiJ6vGc9DE0OWK7Ox6spTPVdQhYch8wGo/kqEadxKg7DDXVtCm4vpkasWn3QS1cRh3pQHGN2nvPR4JU1IdjEtCqRcBw7OTOjHD1elpYPl0rR0GNOvVRKmA3sMl+lWq3m5F2X5Mt93yCCtEHhtowrFUYnJQebq9fEUU0mLQ1mzORkxUU1C4zlFHr+Y6uBG7UIkqoRg1kGTU1Zspjvdc/p2UVjoHSTknbhq5UpyO1b6Bt019OQ9bLoWGMDumDSkHO/C0wORPVdevLCSR1moSMvJCJ+oDfdyB9/8bWaRp1nGfCqEkdZ7xIZRrXQzFqgwzFGhSB1EL2CiOx3HD/pdlT8mrUytAWNtRK+liXlWKoPPnco9HI+dQbgPf9JZuX12YTBkM9OZ6g0yCjLVte6h7q+0eapA+Xs7ck+t9QU416fF1kqE3GRwyYnYm2RfcHhqKbbdLRuAeMRn0AWWZEH8CpWwKfeRI47NPEmchNeCGM2hQqxxpqjUXveVz2OC0DAHZ5TTQi0MvJGEmtHWzrcyTnl5E+DBq1qtMLDshOuhEivQFwKuqDxnjr3w2RBnpd169AMotTGX81+UII4PifAQecDOz3zrasQkMCjRNeuCn0AuzLShnXa/9f9hwutthmqDfRpt0XgVrzZXQVM1mMrDGz5L7IUG+7DzCFmfmYjHgsI18To07am/QZ28zhhFGPtvPUgwqcSylUh9401PqDx0140cGtpzEyne9YohG9Uek0Xf24ni/FtK3MOht7w7WhacKoGemDYtIm6XqmDDUzAjCxA+pA0aFYwokXR4tCmaC3w0mXAJ9flG5v2k55w/MggE89DHz8Xq3enuzNKH3E3weHgRPOy56nyx/cDi9JPh6Tc/SHdXQl2gy5Gf2dvJgFsNdxwLHfBd78o2yepjok5XDGySB9pHaiR1puW/E0k018Ll2YSUpgkxdEfxfdKFZt+LB+efbZoYZ65dPZuqfq6eEApHMLFOiiZonB9tCokx12htuMXFKNeiIa6lOub//tWisDYAz1JnxMZvIWNcgdrm2splo27c2wOhUfrYxl/Dt94GydLjHAg1mjnMrUoCNyLP+QjwNv+5W2hsRk+7CWnX2ltw2VPlzheYTRCBG9APX4U9MiQaa6mZyJQNYnALSlAYCP+kjyMcwS1KHfv/Ur22W3xqNr5aIy0oW06wGY+zu7zKs+YmMYNfd9xZNMPvG5bzkzGskp6Iy6KEY2idpx3fKsRs2t2mhjy+pabbpy4nMiBl8xe9r/bBFPCaNW0sdwFOK63zuBfU4gGvVENNTb7A3s87bob9qQnCebY9SsodZWwVJI/e2YomqLk2QlGkajpp3M2uk0pqzScXU3DfcSx4f2gLz2C8Dux6ZXH9SNw8l/Al79uWwZxjpSRu3SqDXnrU/+1rIN0R0pQ83UQY8bl622quFk0A7H6uiq9jmt8Si/5L4ZfCaZpQMMhpr93SF9JN+16+VilfUZsHqkSxWGujEQPTcso2bYs82hr9oqs0EGk4a2gTLUGenDZqhJ1MfAcOSrevOPojro93RCGmoAyUORMdSMpMFKHxaD7uNM5B7KYcu6zdziUakyVSebEgXqK9gMdeL40KIVuJdJMkORykJx2iaj1+sOEh2zDiax3q44Wosz0Rb1YZtkY8Pbzm+PbHwYNfdQZwx1GelDu8aUoY6nGCvDY2Lk1vA8B9hVF8EYamIQ6chQv079pS5bwHSyJnReiAYwaTPgmXnAvb9PH+MYtdVQezBq1R/pS4AyagU6L0OHMtT7/Uv0ud3+5rQTNjyPrgyn0LRIGgoj0/l0qi/7OBO5h9I11VrHQDztm3vYX/9lrUzbZBON6aiydScpDVeidVZ5c45VnVFTqE6dKsMXWh1su7IkL56cEyp2P6Y9Fd+kUesjH+76dOOtL3OaIdD02h2Oz+ZoO01zLGbUk7RjHKhGbTLUnKHXoj70ulHDTPvtZjuQbLRzU31FRmz8PVcCr2Mckz4QIrofT8/JHuNkKauhbnikMTHquE/o27Pp3zmottj9mGjJXdv08InLqGPQhtx6r2wayiRHprUdABxMEkfqIclpqKlHu6GkD2UsDU5Im6HSnYSKzafW2CWMmrI21S7cbDTVCTljqodNuQypOr7TYcDhn402Ck7Kt7AMU9iaDzIOIe2+vfvyaE0HBe6h1hnl+Hp4R334RKhkNOq4DUzrWFDpoyij1mUFaqhVW+/x5mi0pOvQep1VvRVUnXY4yL6s6MzdzceAtpGk4PqedYSZR6MmbF1F/yimnfgELIbaZwcbhYlrqBWjJh13u/2Ad1+W/o2mGdkkzWBMcbYUrjVAbB5pKrWotXxNizn5QJc09HV9k+OESVNDrRsNU305r7c+JHUZajWj8rB/Bw4/Pd2GtgdK6aY+u4SYwEkfOx2aTsMZ6td9ETj4Y9GwfsNai/ThoVHTa1T3ad2y6PxkEXvD2shcmXseB7zzYj49rU+DM9QG6WO3o4GP3BEZXlMddMJxpMaiTf3ghF9Ek4NsMPl2XKOdDJSfxyONMsh7HR9FFU2No08UCaHTwxX+7Rbg1L+l0/igZkPdmzu86ODeeFNm2NOMTE8zajWlljJQCld4no1RZwz1ANKMuoCh1t/6+ipktI6mmXi2sk26NhDFsvpik+3SO7HosBnqIebF4wu6yL41Goe519O3iYz1kvujmabGqA9du/Vg1DseEi1zAACL741eIEMOQ02lDwB468+yybgQSl1a042KyZmoGKJt5KDyOfyzwJ5vydaTwmeih4lRc8bN1mcU+bIyaiJ9TN48iipSUC8tk/Sx1e7t9YLykIgJa6ht6xXTG0VZIdWoTbofNfguZyLnpVagxrAxiNSa0nmGUQrKKI9M12Z4cfGsZNKEXgeAZwYn/R64+6J2nKsO2zA3D2wPsTIe3DRnJ8gLassXmZPaMDzFzqizonU2D/WAztgtalP1kD/3cLyjjEv6MMhWOj58u2V3IKZ/ZSaFxXX0cZomazB7DrZ9ZizqPg8dqVmskyNZz6Y/qza0pUnmK8SGmjJ0NVrMrG7I1CsPufJtr4LoYenDAtti9UA2PI+u9ayM3itOI/mUcCZy0of+IBVh1KPxFOjJm/PSh4JJ+rB1uK12B17zH/x1ll3aMinf0r2GShhqOpLweUj2PoGpw9ToxUcXZUrK8Yj6UH1xZHqkuU6b2Wb4KY3aJH0wjDpThmXn8e0PaJdvgh6Jos6j+SjQnbdpHhQ+0Q46wVETaGgZKo3VUK93l2lyqCflxO2k+h03YlfEb8au5nI6jN5l1CaWCJjjhRVoeB516E3fGjh9fnZIpjNvjmlaNWoSNphsCOBwJtqwnjPUFkZt2q/RR2s78sv8tN2qsPVewOJ72t+t1+MC479416XA84/yyU3SzPCU6IE1MU2fJVk5h9gm20XxyqLhsX+fpZ8n9TAZagG88VvA7PdkIzk+cifwvf3idIR1Z9YwYaQP330rfaQP3bGn90W9DPXiLsuo9bU4AMY2xM+wmvDFXefgMHDSH4Bt9jGX02H0rqG2kTr6cLhmJk6dEQXcb9DWkuaG9wPD0ZKg+5/IOxytjJox1LqzpxCjVvs/btZmAhwzU5sk0Fldecp+5YfS3/c/Cbjjl/51dWHaVsBitF+yqi03FDDU3IJOO78q+pcHQ9RQ26QOh0atH1MveTHQZtem4b+P9GGUF0TELl9wQPbQFju1/6YvbPpCOvD97b8P/QSw8ingpe/KlsXBNnX6xW+IPnW2r/fF1FIRlogO5TxUfd9nrQ/TyIAuBGUKLNjl1eYyuoDeNdQKXAemb/HtXwbce0n7+/A04OUfAP58evR9+rbAc4+4y2o0gFP/aj5uDc+jhprMfCqiUSu4pI9DPh518Je+m69DHu+1wpu+Dxz7vfznmXDc2cADl7edbYpBjZXRqHOEsnEYnhbvaB/fO5f0wWrUI9ljaqU20YgcchtWt2faUtBFmTiY2K2vRMVtSqFARxtTZ0SRHL5l2Zxo/3JB9Jmasq8zal1qVBOhiKZ86g1tZ+AuRwBX/yewx1vghBot05ExlYiKrrPdYfSwRm3pwPSte9AHgQ/erOmF04CDPgB89hlg9nuBN36zmirZZiYaNWpLiJwvXIZ6aDJw6CfNa30UCg3UJ1NUgKlbAgdoLxIlsxSZoWhbIjUPEofmaj4//fuex/HGatDCqIenRL+/9F1uLTcPo/a97pm7I9nqDTBHfXihgKFWyKxWiGiEse/b27+bGPW2+7SjkLbZO3qxbM+MIChefmokCx3wr+nfM4a6JFflRjM1oPcZNQd1U9WUTiGArfeI1rC99w/th394SrTjODeFughsUR/6VNtNto8YgmiZWe0BJ7eZlwuTNm0/XIOTgJOvBhb83X2e0ti6OYzb7Whe69vpVcBRX00/rN4wxNjnheonfzLsHK3a/LB/jzY45hx+CSPTDXV8X3XHmQmm0Eq2DHXOALz2M/zgTdHnVZ+PPk3Shw9KORM1uUH5aj56V1p+9Jl1mAcDQ8DL3pv9nfbFsgv+v+eqcnMBPNG7htoUyaDwgRuzaw1vu2/0j8I2nz8Ppm1tPvbGb0YG9dazgB1eDqxa3J6dBmQN9bHf9S9X1f8NX4tmAG61u33dAYVt9gJOX1DrpptOvOM3/O9CRKOeIqiMURMpixojZYQmbxHdxybTF9nFu2JD7dPuRaI+GoOx1OYw1HRlPaMW74EXHwXMOjSayKM7hZX0M23rqG9y0AlOoh0b1qWpec0MbLsP8In7gUevBf54GrDVHuXyq8q2ONC7htoF2xrKdcHWwUemR4Z0xouBfd8BXBAv5LLTYcBtZ1dT35efmv+cbhrpItj+QHeapfGGptvtV64sKrvQKKBDPxE5pWfHw2fu/tskC9MGyzoUq/QJr1Ogm1G4QKWPIpi0CXDy5cCdvwH+oL1g1UvkUw+1f3vjN9OjxWFGo86sM8+sAlgXNtkuWqr0Jcf0zfPhZaiFEEcB+C6AAQBnSyn/t9ZapWB5EPLglOurmT20778Ad/06cmw8eWv6mBBtD/quR0Z68h7/BPz7Y+0prHnw0bsiBjNR8Pmlfg+qeshfcky58nRG/f7rsuGJI9OBo76i/cAZRuWgYxaj2nxHdx0OOi16YVAt1QZufQ8buFHdNnsDB/2bf5kK+70jcpB+OR5dcgz4Ze9Lf9cZ9VH/A1zxGfMaKR2QEaLyRN8YacDDUAshBgD8EMDrADwF4DYhxKVSyvvqrVpFky4UfKQCH/zTD6PO5rrJerhbESMNRBvS0k1pN2b47pLxxm9GUT15prpz0A21vpCTCcoYH3By+zf10tBZ437vjF7SsxmNlGJw2DxS2mpPYMm92d/zxMcD7fA2XdP/wI1+53LQw+O84qi1tnnZ+7KGHAB2fxOw8I7iy6oeeCpw65nFzu0D+DDqAwE8IqV8DACEEBcA+CcA9Rrq/U8Cbj8P2Lm34hnRaPTVm7g2vPrzwCPMBqSdwJQtIj9AWej6pCnOWYcQwOcWpY3T9gcCr/pM2nE1MBRFIpXFe64A1jyb/f3lHwCu/4o9XFTHKz8czSF4eUGfAIf93gnceb5fZJDN6fua/4yM877viF6ARSddHf216F9RvPwDbUmtByGkTWMDIIQ4HsBRUsr3xd9PAvByKeWHSLpTAJwCADvssMMBCxYsqKfGAQFVYsVTwMqFwAs9tPFegWmBqE6i1YxCUq2r3Wm49afAdi/1C62boBBCzJVSzuaO+TBqizCn/SDlWQDOAoDZs2dXJCwHBNSMTbe3Lwjfi+i2kQYiltzwNNJAevZjQG74uFifAqDHwW0PYGE91QkICAgIoPAx1LcB2FUIsZMQYhjA2wFcWm+1AgICAgIUnNKHlHJcCPEhAFciCs87V0rJuKIDAgICAuqAVxy1lPJPAP5Uc10CAgICAhj08KJMAQEBAQFAMNQBAQEBPY9gqAMCAgJ6HMFQBwQEBPQ4nDMTC2UqxFIARacmzgDAzJvdqBGueWIgXPPEQNFr3lFKOZM7UIuhLgMhxBzTNMqNFeGaJwbCNU8M1HHNQfoICAgI6HEEQx0QEBDQ4+hFQ31WtyvQBYRrnhgI1zwxUPk195xGHRAQEBCQRi8y6oCAgIAADcFQBwQEBPQ4esZQCyGOEkI8KIR4RAjxmW7XpyoIIc4VQiwRQtyj/baFEOJqIcTD8efm8e9CCPG9uA3mCSFe2r2aF4cQ4oVCiOuEEPcLIe4VQnw0/n2jvW4hxCQhxK1CiLvia/5i/PtOQoh/xNd8YbxUMIQQI/H3R+Ljs7pZ/zIQQgwIIe4QQlwef9+or1kIMV8IcbcQ4k4hxJz4t1r7dk8Yam0D3TcA2APAO4QQe9jP6hv8HMBR5LfPALhGSrkrgGvi70B0/bvG/04B8OMO1bFqjAP4pJRydwAHATgtvp8b83WPAjhCSrkvgP0AHCWEOAjAVwF8O77mZQDU5orvBbBMSvkiAN+O0/UrPgrgfu37RLjmV0sp99Pipevt21LKrv8D8AoAV2rfzwBwRrfrVeH1zQJwj/b9QQDbxn9vC+DB+O8zAbyDS9fP/wD8EdEu9hPiugFMAXA7gJcjmqE2GP+e9HNE67u/Iv57ME4nul33Ate6fWyYjgBwOaKt+zb2a54PYAb5rda+3ROMGsALADypfX8q/m1jxdZSymcAIP7cKv59o2uHeHi7P4B/YCO/7lgCuBPAEgBXA3gUwHIp5XicRL+u5Jrj4ysAbNnZGleC7wD4NIBW/H1LbPzXLAFcJYSYG2/qDdTct702DugAvDbQnQDYqNpBCDENwMUAPialXCnMm7JuFNctpWwC2E8IsRmA3wPYnUsWf/b9NQshjgGwREo5VwhxuPqZSbrRXHOMg6WUC4UQWwG4WgjxgCVtJdfcK4x6om2gu1gIsS0AxJ9L4t83mnYQQgwhMtLnSykviX/e6K8bAKSUywFcj0if30wIoQiRfl3JNcfHNwXwfGdrWhoHA3iTEGI+gAsQyR/fwcZ9zZBSLow/lyB6IR+Imvt2rxjqibaB7qUA3h3//W5EGq76/V2xp/ggACvUcKqfICLqfA6A+6WU39IObbTXLYSYGTNpCCEmA3gtIgfbdQCOj5PRa1ZtcTyAa2UsYvYLpJRnSCm3l1LOQvTMXiulfCc24msWQkwVQkxXfwM4EsA9qLtvd1uY10T2owE8hEjX+1y361Phdf0GwDMAxhC9Xd+LSJe7BsDD8ecWcVqBKPrlUQB3A5jd7foXvOZDEA3v5gG4M/539MZ83QD2AXBHfM33APjP+PedAdwK4BEAvwMwEv8+Kf7+SHx8525fQ8nrPxzA5Rv7NcfXdlf8715lq+ru22EKeUBAQECPo1ekj4CAgIAAA4KhDggICOhxBEMdEBAQ0OMIhjogICCgxxEMdUBAQECPIxjqgICAgB5HMNTjzt1VAAAADElEQVQBAQEBPY7/D/HEAGSq+AvcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f_max = abs(v_max.dot(M).dot(v_max))\n",
    "print('f from the product of M and v_max = ' + str(f_max))\n",
    "print('The largest (absolute value) eigenvalue = ' + str(lm_max))\n",
    "\n",
    "num_v = 500\n",
    "#  随机建立多个归一化向量\n",
    "vecs = np.random.randn(num_v, dim)\n",
    "vecs = np.einsum('na,n->na', vecs, 1/np.linalg.norm(vecs, axis=1))\n",
    "# 计算每个向量的f\n",
    "f = abs(np.einsum('na,ab,nb->n', vecs, M, vecs.conj()))\n",
    "\n",
    "# 画图展示由最大本征值给出的\n",
    "x = np.arange(num_v)\n",
    "y = np.ones(num_v, ) * f_max\n",
    "plt.plot(x, y, '--')\n",
    "plt.plot(x, f)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The dominant eigenvalue and eigenvector by eigs:\n",
      "[-4.83004377+0.j]\n",
      "[ 0.13390826+0.j -0.32371695+0.j  0.17090423+0.j  0.92090589+0.j]\n",
      "\n",
      "The dominant eigenvalue and eigenvector by eig:\n",
      "-4.830043768468989\n",
      "[-0.13390826  0.32371695 -0.17090423 -0.92090589]\n"
     ]
    }
   ],
   "source": [
    "# 使用scipy中的eigs仅求最大几个本征值与本征向量，节省计算量\n",
    "import scipy.sparse.linalg as la\n",
    "\n",
    "lm1, v1 = la.eigs(M, k=1, which='LM')\n",
    "print('The dominant eigenvalue and eigenvector by eigs:')\n",
    "print(lm1)\n",
    "print(v1.reshape(-1, ))\n",
    "\n",
    "print('\\nThe dominant eigenvalue and eigenvector by eig:')\n",
    "print(lm_max)\n",
    "print(v_max.reshape(-1, ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eig0(mat, it_time=100, tol=1e-15):\n",
    "    \"\"\"\n",
    "    :param mat: 输入矩阵（实对称阵）\n",
    "    :param it_time: 最大迭代步数\n",
    "    :param tol: 收敛阈值\n",
    "    :return lm: （绝对值）最大本征值\n",
    "    :return v1: 最大本征向量\n",
    "    \"\"\"\n",
    "    # 初始化向量\n",
    "    v1 = np.random.randn(mat.shape[0],)\n",
    "    v0 = copy.deepcopy(v1)\n",
    "    lm = 1\n",
    "    for n in range(it_time):  # 开始循环迭代\n",
    "        v1 = mat.dot(v0)  # 计算v1 = M V0\n",
    "        lm = np.linalg.norm(v1)  # 求本征值\n",
    "        v1 /= lm  # 归一化v1\n",
    "        # 判断收敛\n",
    "        conv = np.linalg.norm(v1 - v0)\n",
    "        if conv < tol:\n",
    "            break\n",
    "        else:\n",
    "            v0 = copy.deepcopy(v1)\n",
    "    return lm, v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The dominant eigenvalue and eigenvector by the iterative method:\n",
      "4.830043768468986\n",
      "[ 0.13390826 -0.32371695  0.17090423  0.92090589]\n"
     ]
    }
   ],
   "source": [
    "lm2, v2 = eig0(M)\n",
    "\n",
    "print('\\nThe dominant eigenvalue and eigenvector by the iterative method:')\n",
    "print(lm2)\n",
    "print(v2.reshape(-1, ))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
